; docformat = 'rst'
;
;+
;
; :Purpose:
;	  Format measurements for use in CELS framework. 
;
; :Inputs:
;	  Sim_name: String variable containing the name for the simulation
;
; :Requires:
;   measurements.inp.sav file.  File must contain:
;   
;     *mf: mole fractions with dimensions of [station, time]. For missing data use !values.f_nan.*
;   
;     *mf_error: repeatability of each measurement. Dimensions of [station, time].*
;   
;     *mf_time: Julian time for mid-point of each observation averaging period, with same interval as LPDM. Must span entire date range.*
;   
;     *mf_lon: longitudes of stations*
;   
;     *mf_lat: latitudes of stations*
;   
;     *mf_station_ID: station codes, with same code as in param file*
;
; :Keywords:
;	  Afternoon:  Set this keyword to use only observations between 1000 and 1600hrs local time
;   Scale:  Fractional change to the measurement scale.  E.g. to increase by 2% use Scale=0.02
;   Spinup: Process measurements_spinup.inp.sav
;   
; :Outputs: 
;   Creates measurements.sav file in Data/{sim_name}/ directory or 
;   measurements_spinup.sav if Spinup keyword set
;
; :Example::
;	  IDL> cels_measurements, 'SF6', /afternoon, scale=-0.02
;
; :History:
; 	Written by: Matt Rigby, MIT, May 29th 2011
;
;-
Pro cels_measurements, sim_name, afternoon=afternoon, scale=scale, spinup=spinup

  compile_opt idl2
  on_error, 2

	;Set parameters for pollutant
	;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  if keyword_set(spinup) then spinup_str='_spinup' else spinup_str=''
  print, strcompress("CELS_MEASUREMENTS: formatting measurements" + spinup_str)
  
  StartY=cels_get_parameter(sim_name, 'StartY')
  EndY=cels_get_parameter(sim_name, 'EndY')
  if keyword_set(spinup) then startY_IC=cels_get_parameter(sim_name, 'StartY_IC')
  Stations_all=[cels_get_parameter(sim_name, 'Stations'), cels_get_parameter(sim_name, 'BG_Stations')]
  nStations_all=n_elements(stations_all)
  
  ;Get measurement input file
  ;MUST CONTAIN:
  ; mf: mole fractions with dimensions of [station, time]. For missing data use !values.f_nan.
  ; mf_error: repeatability of each measurement. Dimensions of [station, time].
  ; mf_time: Julian time, with same interval as LPDM. Must span entire date range.
  ; mf_lon: longitudes of stations
  ; mf_lat: latitudes of stations
  ; mf_station: station codes, with same code as in param file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  restore, cels_filestr(/Input, sim_name + '/measurements/measurements' + spinup_str + '.inp.sav')
  
  ;Select only required data range
  if keyword_set(spinup) then begin
    wh=where( (mf_time gt julday(1, 1, StartY_IC, 0)) and (mf_time lt julday(1, 1, StartY, 0)), count )  
  endif else begin
    wh=where( (mf_time gt julday(1, 1, StartY, 0)) and (mf_time lt julday(1, 1, EndY, 0)), count )
  endelse
  
  if count eq 0 then begin
    message, 'Invalid date range'
  endif else begin
    mf_time=mf_time[wh]
    mf=mf[*, wh]
    mf_error=mf_error[*, wh]
  endelse

  TimeSize=n_elements(mf_time)

  ;Sort measurement stations
  mf_sorted=fltarr(nStations_all, TimeSize)
  mf_error_sorted=fltarr(nStations_all, TimeSize)
  mf_lon_sorted=fltarr(nStations_all)
  mf_lat_sorted=fltarr(nStations_all)
  for si=0, nStations_all-1 do begin
    wh=where(mf_station_ID eq stations_all[si], count)
    if count eq 0 then begin
      message, 'Missing station ', stations_all[si]
    endif else begin
      mf_sorted[si, *]=mf[wh[0], *]
      mf_error_sorted[si, *]=mf_error[wh[0], *]
      mf_lon_sorted[si, *]=mf_lon[wh[0]]
      mf_lat_sorted[si, *]=mf_lat[wh[0]]
    endelse
  endfor
  mf=temporary(mf_sorted)
  mf_error=temporary(mf_error_sorted)
  mf_lon=temporary(mf_lon_sorted)
  mf_lat=temporary(mf_lat_sorted)
  mf_station_ID=stations_all

  ;Specify where finite measurement locations are (time and station)
  ;If afternoon keyword is set, only keep afternoon measurements
  if keyword_set(afternoon) and (keyword_set(spinup) eq 0) then begin
    local_hour=fltarr(nStations_all, TimeSize)
    dummy=0
    caldat, timeconc, dummy, dummy1, dummy2, hour
    for si=0, nStations_all-1 do local_hour[si, *]=fix(float(hour) + 24.*mf_lon[si]/360.) mod 24
    mf_wh=where(local_hour gt 10. and local_hour lt 16. and finite(mf), count)
  endif else begin
    mf_wh=where(finite(mf), count)
  endelse
  
  if count eq 0 then begin
    message, 'No data in range'
  endif else begin
    mf_ai=array_indices(mf, mf_wh)
    mf=mf[mf_wh]
    mf_error=mf_error[mf_wh]
  endelse

	;Scale observations if required (useful for estimating error due to calibration scale)
  if keyword_Set(scale) then begin
    mf=mf*(1.+scale)
  endif

  save, filename=cels_filestr(/Input, sim_name + '/measurements/measurements' + spinup_str + '.sav'), $
    mf_time, mf, mf_lon, mf_lat, mf_wh, mf_ai, mf_station_ID, mf_error

end